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ABSTRACT 

Aims. We explore the effects of an outer stably stratified coronal envelope on rotating turbulent convection, differential rotation, and 
large-scale dynamo action in spherical wedge models of the Sun. 

Methods. We solve the compressible magnetohydrodynamic equations in a two-layer model with unstable stratification below the 
surface, representing the convection zone, and a stably stratified coronal envelope above. The interface represents a free surface. We 
compare our model to models that have no coronal envelope. 

Results. The presence of a coronal envelope is found to modify the Reynolds stress and the A effect resulting in a weaker and non- 
cylindrical differential rotation. This is related to the reduced latitudinal temperature variations that are caused by and dependent 
on the angular velocity. Some simulations develop a near-surface shear layer that we can relate to a sign change in the meridional 
Reynolds stress term in the thermal wind balance equation. Furthermore, the presence of a free surface changes the magnetic field 
evolution since the toroidal field is concentrated closer to the surface. In all simulations, however, the migration direction of the mean 
magnetic field can be explained by the Parker-Yoshimura rule, which is consistent with earlier findings. 

Conclusions. A realistic treatment of the upper boundary in spherical dynamo simulations is crucial for the dynamics of the flow and 
magnetic field evolution. 
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1. Introduction 

The Sun has an activity cycle of about 11 years, with an 
underlying magnetic field that oscillates with a period of 
around 22 years. A dynamo operating in the convection zone 
below the solar surface i s responsible for generating cyclic 
magnetic fields (se e e.g., Brandenburg & SubramaniarJ 120051 : 
Charbonneaul 120051 and reference therein). The occurrence of 
sunspots, which is the main surface manifestation of the solar cy¬ 
cle, varies regularly over the cycle. At the beginning of each cy¬ 
cle sunspots tend to appear at mid-latitudes, while toward the end 
of each cycle they tend to appear at low latitudes. It is believed 
that these sunspots and their occurrence are connected to an un¬ 
derlying toroidal magnetic field, which is migrating equatorward 
during the cycle. Theoretical models of solar magnetic field evo¬ 
lution have been studied for several decades. Mean-field models, 
where turbulence e ffects are parameterized through transport co¬ 
efficients (see e.g.. lKrause & Radlerlll98Cih . have been success- 
ful in producing some observed magnetic field prope rties (e.g., 
iKanvlh et al.l 120061 : iKitchatinov & OlemskovI l2012h . Another 
class of dynarno models relies on the Babcock-Leighton effect 
(lBabcockll96ll l Lei ghtoni 19641) and flux transport by meridional 


tion, where the equations of magnetohydrodynamics are solved 
directly. However, the increasing computing power has led to the 
successful reproduction of some observed features of the solar 
magnetic field by such models. 

For a long time global simulation s were only ab l e to gener¬ 
ate ei t her poleward migra ti ng fields ( Gilm^ 19831: iBrun et ahl 
12004 iKanvla et all 120101: iNelson et al.1 12013h~ or oscillatory 
ones 


IKanvla et al.l IJNelson et al.l I/BI JI) or oscillatory 

with no cle ar migration pat tern (i Ghizaru et alJ 20Tft 


Racine et al 


Kanvla et al 


20nk For the first time. lSchrinner et al.l (1201 1 ) and 


(12012l) could produce clear equatorward migration 
of the toroidal magnetic field. Stra tification and rotati on rate had 
to be high enough for this to work (iKanvla et al.l2013h . Recently, 
several groups have been able to pro d uce grand minima-type 
events (IPassos & Charbonneau 201^ Augustson et alJ 120151 : 


circul ation (e.g., Choudhuri et al.in~995l:lDikpati & Charbonneaul 

Il999l) . The low Reynolds numbers compared with the Sun limit 
the usefulness of global simulations of self-consistent convec- 


iKapvla et alJ 2016a ). In Kapvla et al. (l2016a ). a secondary dy- 
namo mode dis t urbed the surface field of the primary mode. 
IWarnecke et al.l (l2014t) co uld explain the equatorw ard migra- 
tio n seen in simulat i ons of iKapvla et al.l (12012i 12013l) and those 
of lAugustson et all (12015l) as a propa gating afl dynamo wave 
follow ing the Parker-Yoshimura rule (lParkeil[T955l : lYoshimui^ 
119751) . Here, a was estimated via the kinetic helicity and £2 is 
the local solar rotation rate. An equatorward migrating dynamo 
wave is possible i f a is positive (negative ) in the northern (south¬ 
ern) hemisphere (ISteenbeck et al.l[l966l) and the radial gradient 
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of Q is negativ e. These interpretatio ns have been confirmed in¬ 
dependently by iDuarte et al.l (12016h and by the computation of 
turbule nt transport coe f ficients obtained through the test-field 
method IWarnecke et al.l (12016bll . 

Given that this simple relation can describe the behavior of 
dynamos driven by self-consistent turbulent convection in sim¬ 
ulations, the Parker-Yoshimura rule can also be a possibility to 
explaining the equatorward migration of the magn etic field of 
the Su n. Indeed, the differential rotation of the Sun (ISchou et al.l 
Il998ll has a negative radial gradient in the near- s urface shear 
layer (iThomnson et al.l ri996ri^rekat et akl l2014l | 2016ll. This 
make s it a possible location of the solar dynamo (iBrandenburgI 
12001 . 


The generation of the solar differential rotation can well 
be described by mean-field models, where the productive 
parts of the o ff’-diagonal Rey nolds stress are parameterized by 
the A effect (iRudigeii Il989l) and the turbulent heat transport 
in terms of anisotropic turbulent heat conductivity (see also 
[Brandenburg et al]ll992h . These models can reproduce the sur¬ 
face differential rotation, the sp oke-like rotation profile and 
the near-surface shear lay er fe.g.. lKitchatinov & Rudig^ll995l 
l2005tlRfidiger et al.l2013ll . With global models of turbulent con- 
vection, it is challenging to generate such rotation profiles (e.g. , 
iDe Rosa et al.l 12001 IBrandenburgIl2007h . iMiesch et al.l (12001 
were able to produce a spoke-like rotation profile by impos- 
ing a latitudinal e ntropy gradient at the bottom boundary, while 
iBrun et a kl (1201 ll) used a s t ably st ratified layer below the convec¬ 
tion zone. iGuerrero et al.l (12013l) could produce a near-su rface 
shear layer at lower latitudes. RecentlyjHotta_etj^ 201 5l) used 
a reduced sound speed technique (iHotta et al. 12012 ) achieving 
high stratification to produce a near-surface shear layer, the gen¬ 
eration of which was suggested to be related with the meridional 
Reynolds stress component Qrg. However, they still had cylin¬ 
drical rotation contours within the convection zone. 


IWarnecke et al.l (l2013al) applied a different approach and 
used an outer coronal envelop e above the dynamo domain 
(IWarnecke & Brandenburgl2010ll to reproduce spoke-like differ¬ 
ential rotation at low latitudes with a weak near-surface shear 
layer. This two-layer approach has also been used to success¬ 
fully simulate coro nal ejections driven by dynam os arising from 
forced turbulence (IWarnecke et al. 201 iL ^ 12 al) as well as by 
convective dynamos (IWarnecke et al.l 2012 bh . l^rthermore. the 
outcome of dynamo simulations suggests that the presence of a 
coronal e nvelope supports the dynamo and leads to a higher field 
strength dWarnecke & Brandenburgl2014li . 


To inyestig a te and follow up on the findings of 
IWarnecke et al.l (l2013ah . we perform a detailed study of 
similar simulations with and without a coronal enyelope to 
inyestigate the effect of a coronal enyelope as a free boundary 
on a conyectiyely driyen dynamo. We yary the size of the 
enyelope, as well as the cooling profile, the magnetic boundary 
condition, and the rotation rate. We analyze the effect on the 
flows, differential rotation, and the magnetic field eyolution. 
Eyen though the solar corona most likely has limited influence 
on the dynamics of subsurface flows and the eyolution of the 
magnetic fields in the Sun, these studies are important for inyes- 
tigating different influences and effects on conyectiye dynamo 
simulations. Eyery simulation, in which we better understand 
the mechanism causing flow and magnetic field eyolution, will 
bring us a step closer toward understanding the dynamics of the 
interior of the Sun and other stars. 
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Fig. 1. Radial profiles of azimuthally and latitudinally averaged tem¬ 
perature {T)g^ (a) and density {p)b^ (b) normalized by their values at the 
bottom of the domain, Tq or po, respectively for Set A. The inlays show 
the radial profiles near the surface. 


2. Model and setup 


Our s et up is similar to the one-layer model ofjKag^laet^ 


( 2 () 1 _ 2 |. l2013h and the two-layer model of IWarnecke et al 
( 2013ah. both of wh ich have recently also been used 


IWarnecke et al.l (I20l4) . We use a wedge in spherical polar coor¬ 
dinates (r, G, (p), in which the layer below the surface (ro < r < R) 
represents the convection zone. Here, R is the solar radius and ro 
corresponds to the bottom of the convection zone at r = 0.7 R. 
The layer above the surface (R < r < Rc) represents a sim¬ 
plified coronal envelope, which extends to outer radius Rc- The 
domain spans 15° < 6 < 165° in colatitude and 0° < <p < 90° in 
longitude. We solve the equations of compressible magnetohy¬ 
drodynamics. 


dA 

— ^uxB- Pc'?/, 
ot 


( 1 ) 


Dlnp 

Df 


-V M, 


( 2 ) 


Dm 1 

-— - g - 2r2o xu + -{JxB - Vp -i- V ■ 2vpS), (3) 

Df p 

= --V ■ + 2yS2 + ^72 _ ( 4 ) 

Df p ' ' p 

where the magnetic field is defined via the vector potential B - 
Vx A, ensuring the solenoidality of B at all times, J - VxB/po 
is the current density with po being the vacuum permeability, 77 
is the magnetic diffusivity, v is the kinematic viscosity, u is the 
plasma velocity, p is the mass density, s is the specific entropy. 
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and D/Df = djdt + « • V is the Lagrangian derivative. The trace¬ 
less rate-of-strain tensor is given by 

S/y — "b 

where sem i colon s denote covariant differentiation; see 
iMitra et ^ (l2009l) for details. The gravitational accelera¬ 
tion is given by 

g - -GMrjr^, (6) 

where G is Newton’s gravitational constant and M is the mass 
of the Sun. In addition Jlo = f2o(cos 0, - sin 6, 0) is the rotation 
vector, where Qq is the rotation rate of the comoving frame. Us¬ 
ing the ideal gas law, the pressure is given hy p - {y - l)pe, 
where y = cp/cy = 5/3 is the ratio of specific heats at constant 
pressure and constant volume, respectively, and e - cyT is the 
internal energy density, which is related to the temperature T. 
The two diffusive heat fluxes are defined as 

/rrad ^ = -XsGspTVs, (7) 

where is the radiative heat flux with the radiative heat 
conductivity K and is the subgrid scale (SGS) heat flux 
that carries the unresolved turbulent heat flux of conve ction 
with the SGS heat diffusivity ;j'sgs ; see lKanvla et al.l (12013h and 
IWarnecke et al.l (12013ah for details. Finally, the function Tcooi re¬ 
laxes the temperature toward a predefined profile T^efir) 


Tcooi — ro/(r) ^ r ’ (^) 

Tref(r) 

where Fq is a cooling luminosity, /(r) is a profile function tend¬ 
ing to_unit^Mri_?^ R and going smoothly to zero in r < /?, 
see IWarnecke et all (l2013ah for details. Figure [1] shows the cor¬ 
responding temperature and density stratifications for the runs in 
Set A, see Table[T] 

We use is entropic, hydrostatic initial conditions, as i n pre¬ 
vious models dKanvla et al Jl20 1 ^ IWarnecke et al.llT013ah . This 
initial setup is not in thermal equilibrium, but the flux at the 
lower boundary exceeds the flux leaving at the outer boundary, 
resulting the onset of convective instability. Furthermore, we ini¬ 
tialize the magnetic field as a white noise seed field in the con¬ 
vection zone. We apply periodic boundary conditions in the az¬ 
imuthal {(p) direction. For the velocity field we apply stress-free 
boundary conditions at the radial and latitudinal boundaries. The 
magnetic field follows a perfect conductor condition at the lower 
radial and at the latitudinal boundaries. We force the field to be 
radial at the top boundary. Furthermore, the temperature gradi¬ 
ent at the bottom boundary is fixed to have a constant heat flux 
into the domain, and the latitudinal boundaries are impermeable 
for heat fluxes. On the upper radial boundary we either apply a 
black body condition, 

= -K^ -XSGSPT^, (9) 

or or 

or a constant temperature 

T = = Rc). (10) 


In the former case the heat is transported out o f the domain via 
an enhanced SGS flux; see 
case via a cooling flux; see 


Kag^lajet^lJj 20131) and in the latter 


Warnecke et al.l ( 2013ah . 


We apply roughly 10-50 times larger values of viscosity and 
magnetic diffusivity in the coronal envelope compared to the 
convection zone to avoid high velocity amplitudes and strong 


shear flows resulting in small-scale magnetic field enhance¬ 
ments. For the transition region we use a hyperbolic tangent 
function with a radial width of 0.01 R 2 A r - 1.06 R, which is 
the largest value feasible. In the case of Run A5, the viscosity in 
the coronal envelope is just three times higher than in the convec¬ 
tion zone, but we apply an additiona l sho ck viscosity and sh ock 
diffusivity; see iHaugen et al.l (l2004l) and iGent et akl (1201 3l) for 
details regarding their implementation. The heat conductivity K 
in the coronal envelope is chosen such that the heat diffusivity 
X = KIpcp is constant. 

We characterize the runs by the values of the input parame¬ 
ters: Prandtl number Pr = vfx'", sub-grid scale Prandtl number 
PrsGS = ^/Tsgs’ nt^gnstic Prandtl number PrM = vlg, Taylor 
number Ta = {2Q.q{R - ro)^/v)^, and Rayleigh number. 


GM{R - rp)^ / 1 d^\ 
^4 gs^^ I cp dr /,.,0 85 s’ 


( 11 ) 


which is obtained from a hydrostatic one-dimensional model for 
the same initial setup with x"' - xif - 0.85 R) and x'sgs ~ 
XSGsif - 0.85 R). Furthermore we define the fluid and magnetic 
Reynolds numbers. Re = Utmslykf and ReM = Ui-mslgkf, respec¬ 
tively, the Coriolis number Co = 2Qo/Mrms^f and the Peclet num¬ 
ber Pe = Mi-ms/TsGS^f’ where k{ - 2 jtI{R - ro) » 21/R is used as 
a reference wavenumber and where Urms is the typical turbulent 
velocity in the convection zone defined as 


ns — '^3/2 (m/ -(- ll‘Q)0(pr<R i 


( 12 ) 


which corrects for the removal of the differential rotation- 
dominated 1 / 0 . The slow mean meridional flows, however, are 
not removed. Azimuthal averages combined with time averages 
in the saturated stage are referred to as mean and indicated with 
an overbar, e.g., B, while other averages are indicated as {.) with 
the spatial directions as indices. The index 0 refers to the value 
at the bottom of the domain, that is po and Tq. We also use the 
meridional distribution of turbulent velocities 

(13) 


where the fluctuating velocity is defined via u' — u -u. Thus, 
meridional mean flows are here removed. The mesh is equidis¬ 
tant in all directions. The grid resolutions are given in Table[T] 

We express o ur _ results in physical units following 

iKapvla et al.l(l2013ll2014l) and IWarnecke et akl (12014l) by choos¬ 
ing a normalized rotation rate Q =Qo/Go, where Qq = 2.7 x 
10“^ s ' is the solar rotation rate. The simulations were per¬ 
formed with the Pencil Cod^ which uses a high-order finite dif¬ 
ference method for solving the compressible equations of mag¬ 
netohydrodynamics. 


3. Results 

In this work we compare and analyze 13 runs divided into two 
sets based on the rotation rate. In Set A the normalized rotation 
rate is Q = 5 whereas the runs in Set B have a rotation rate 
of Q = 3. For both sets we investigate the effects of a cool¬ 
ing layer and the blackbody boundary condition as well as the 
size of the coronal envelope. A summary of the runs can be 
found in Table [T] and the stratifications of temperature, density 
and entropy of Runs A are sho wn in Figs. [J and |3 Run Al is 
nearly the same as Run B4m of iKapvla et air(l2012l) . Run Cl of 

' http://github.com/pencil-code/ 


Article number, page 3 ofll7l 



















































A&A proofs: manuscript no. paper 


Table 1. Summary of runs. 


Run 

Resolution 

Q 

Rc/R 

Pr 

PfSGS 

PrM 

Ta 

Ra 

Po/Psurf 

PO / Ptop 

Co 

Re 

Pe 

Al 

180 

X 

256 

X 

128 

5 

1.0 

73 

2 

1 

1.25 

HF” 

4.0 

4^ 

31 

31 

8.1 

35 

69 

Ale 

128 

X 

256 

X 

128 

5 

1.01 

78 

2 

1 

1.25 

10 * 

3.9 

10 ’ 

51 

80 

9.4 

30 

60 

Alc2 

128 

X 

256 

X 

128 

5 

1.01 

68 

2 

1 

1.25 

10 * 

2.9 

10 ’ 

22 

31 

9.4 

30 

60 

Alpc 

180 

X 

256 

X 

128 

5 

1.0 

73 

2 

1 

1.25 

10 * 

4.0 

10 ’ 

35 

35 

8.3 

34 

68 

A2 

290 

X 

256 

X 

128 

5 

1.2 

79 

2 

1 

1.25 


4.0 

10 ^ 

59 

7970 

10.4 

27 

54 

A3 

400 

X 

256 

X 

128 

5 

1.4 

79 

2 

1 

1.25 

10 * 

3.9 

10 ’ 

61 

10510 

10.1 

27 

53 

A3t 

400 

X 

256 

X 

128 

5 

1.4 

67 

2 

1 

1.25 

10 * 

3.0 

10 ’ 

23 

5380 

10.2 

28 

55 

A4 

520 

X 

256 

X 

128 

5 

1.6 

79 

2 

1 

1.25 

10 * 

3.9 

10 ’ 

59 

12740 

10.8 

26 

52 

A5 

600 

X 

512 

X 

256 

5 

1.4 

35 

0.5 

0.5 

5.00 


1.8 

10 ^ 

22 

4937 

8.7 

65 

32 

B1 

180 

X 

256 

X 

128 

3 

1.0 

73 

2 

1 

0.45 

10 * 

4.0 

10 ^ 

36 

36 

4.3 

39 

79 

Blc 

128 

X 

256 

X 

128 

3 

1.01 

78 

2 

1 

0.45 

10 * 

3.9 

10 ’ 

52 

83 

5.6 

33 

67 

B3 

400 

X 

256 

X 

128 

3 

1.4 

78 

2 

1 

0.45 

10 * 

3.9 

10 ’ 

58 

8803 

5.5 

31 

61 

B3t 

400 

X 

256 

X 

128 

3 

1.4 

79 

2 

1 

0.45 

10 * 

3.0 

10 ’ 

24 

5380 

5.4 

32 

64 


Notes. The second to ninth columns show quantities that are input parameters to the models, whereas the quantities in the last five columns are 
results of the simulations computed from the saturated state. All quantities are volume averaged over the convection zone r < R, unless explicitly 
stated otherwise, n = is the normalized rotation rate and Rc/R is the outer radius of the domain. Here, pg, p^urt, and ptop are the latitudinal 

and azimuthally averaged density at the bottom (r = 0.7 R), the surface (r = R), and the top (r = Re) of the domain. 


( 201.3h. Run I of lWarnecke et al.l (l2014l) . and the run 


of Warnecke et alJ (1201 6bh . where a blackbody boundary condi¬ 
tion is used. However, we choose a slightly higher stratification 
and a slightly lower value of Prscs^ na mely 2 instead of 2.5 ; 
therefore Run A1 is the same as Run D3 of lKanvla et al.l (l2016cll . 
In Runs Ale and Alc2, the blackbody boundary condition is re¬ 
placed by a shallow {R < r < Rc - 1.01 R) cooling layer, where 
in the former case the temperature minimum is below the surface 
(r = R) and in the lat ter above the surface. T hese two runs have 
been recently used in IWarnecke et al.l (l2014t) as Runs III and IV, 
respectively. The only difference between Runs Alpc and A1 is 
the use of a perfect conductor condition instead of a radial field 
condition for the magnetic field at the top boundary. 


The other runs of Set A have a coronal envelope with differ¬ 
ent outer radii Rc- Runs A2, A3, and A4 have the same cooling 
function as Run Ale, where the temperature reaches a minimum 
below the surface. The temperature increases to a constant coro¬ 
nal value, which is more than twice the value at the bottom of 
the convection zone; see Fig.lTfa). This results in a positive en¬ 
tropy gradient above r - 0.97 R, where the convection ceases 
and drops by a factor of two; see Fig. |2a,b). In Run A3t, 
the same cooling function is applied as in Run Alc2, leading to 
a temperature minimum above the surface at r a; 1.01 R. This 
causes the entropy gradient to become positive at the surface 
(r = R) and an increase of all the way to the surface; see 
Fig. Ha ,b). Already here, we can state that the use of cooling 
profiles in Runs Ale and A3t reproduce most properties of the 
density, temperature and entropy stratification as the blackbody 
boundary conditions in Run Al. In Run A5, we lower Prscs and 
PrM to 0.5 and therefore increase the fluid Reynolds number. As 
the heat flux at the boundary is the same as in the other runs, 
increases only slightly, see Fig. Ha), leading to a reduced 
Coriolis number; otherwise Run A5 is similar to Run A3t. How¬ 
ever, the use of a lower value of viscosity in combination with 
a shock viscosity allows higher velocities in the coronal enve¬ 
lope as shown in Fig. Ha). The runs of Set B are essentially the 
same as the corresponding runs in Set A with a lower rotation 
rate. The radial temperature, density, entropy and velocity pro¬ 


files behave similarly as in Set A and are therefore not shown 
here. In the following we investigate the influence of the coronal 
layer on mass flux and temperature distribution (Section ITTI) as 
well as on differential rotation and meridional circulation (Sec¬ 
tion [321). Furthermore we discuss Reynolds stresses and the A 
effect (Section 13.31) and their contribution to the thermal wind 
balance (Section [3.41) . Then we investigate the influence of the 
magnetic top boundary on the field structure near the surface 
(Section [331 ) and the magnetic field evolution (Section iTbl ). 


3.1. Mass flux and temperature distribution at the boundary 

We begin by inspecting the influence of the top boundary on the 
mass flux and the temperature distribution; see Fig. He). Except 
for Run Al, the mass flux at r = R is nonvanishing, showing 
a strong latitudinal dependency. Near the equator it is positive 
(outflow) and at latitudes around +30° negative (inflow) suggest¬ 
ing a circulation in the coronal envelope. At mid-latitudes the 
mass flux becomes positive again, but is fluctuating around zero 
toward higher latitudes. The flow structure is strongly influenced 
by the rotation, as seen from the alignment with the inner tan¬ 
gent cylinder; see Fig. He). Although the mass flux in the runs 
with an extension above the surface has non-zero values, they 
are small in comparison to pwrms at the surface. Furthermore, 
there is no qualitative difference between runs with a cooling 
layer Rc = 1.01 R (Runs Ale and Alc2) and a coronal envelope 
Rc ^ 1.2R (Runs A2, A3, A3t, and A4). In fact, Run Alc2 has 
the highest mass flux through the surface; see Fig. He). Further¬ 
more, Run A5 is similar to the other runs showing more fluc¬ 
tuations as a function of latitude but with a similar magnitude 
of variations as in the other runs. This shows that the influence 
of a coronal envelope via a radial mass flux is small and can be 
neglected. Furthermore, we find no indication that the viscosity 
profile has a major influence on the mass flux. 

As a second step we investigate the latitudinal temperature 
variation, AT - {T — {T}g)l{T)g, at two radii; see Fig. [3] In 
general the surface perturbations are strongest near the poles, 
decrease toward a minimum at mid-latitudes, and increase again 
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Fig. 2. Radial profiles of azimuthally and latitudinally averaged tur¬ 
bulent rms velocity (a) and {s}g^ (b) in m/s or normalized by cp 
respectively for Set A. Mean radial mass flux p Ur/Psurf“rms (c) through 
the surface (r = R) for Set A. The horizontal line indicates the zero 
value and the three vertical thin lines indicate the equator (d = 90°) and 
the intersection with the inner tangent cylinder (6-90° » ±45° latitude). 


below ±20° latitude. However, in the middle of the convection 
zone the temperature minimum is at the equator. This is a clear 
indication of a strong rotational influence on the temperature dis¬ 
tribution. Run Al shows the largest relative temperature pertur¬ 
bations; over 0.3 near the poles and up to -0.1 at ±20° latitude. 
Here, the radiative boundary condition lets the temperature at the 
surface evolve more freely, which leads to this strong variation. 

Runs Ale and Alc2 have a similar distribution, but with 
around three times smaller values. There the cooling layers cool 
the temperature to a certain latitudinally independent value with 
a relaxation time equal to the turnover time. This leads to a re¬ 
duction of the temperature perturbations near the surface. The 
temperature difference AT is significantly reduced for all runs 
with a coronal envelope. However, the temperature is still higher 
near the poles and lower near the equator than the latitudinal av¬ 
erage. In these runs the coronal envelope with its mass and heat 
capacity serves as a buffer in smoothing the temperature at the 
surface. The influence of the cooling layer and the coronal en¬ 
velope seems to penetrate also deeper down and influences the 
temperature variations in the middle of convection zone. The dif¬ 
ference in the temperature profiles caused by the coronal enve¬ 
lope can also influence differential rotation; see Section [372l 



30 60 90 120 150 

B [deg] 


Fig. 3. Latitudinal dependence of AT = (T - {T)g)l(T)e at radius r = 
0.98 R (a) and at radius r = 0.85 R (b) for Set A. Linestyles as in Fig.[T] 
The thin black lines indicate the zero value, the equator (6 = n/2) and 
the location of the inner tangent cylinder. 


In the higher Re and thus more turbulent Run A5, the temper¬ 
ature variation is similar to that of the other runs with a coronal 
envelope. The runs with slower rotation (Set B, not shown here) 
show a similar behavior, but the latitudinal temperature pertur¬ 
bations are weaker than in the more rapidly rotating runs, due to 
the reduction of the rotational influence. 

On the solar surface, systematic latitudinal temperature vari¬ 
ations have not been observed. However, variations in the range 
of a few kelvin, so less than 0.1% of the surface tempera¬ 
ture, are below the measurable range. This value is exceeded in 
all of our simulations. The temperature difference between the 
pole s and equator can influence the differe ntial rotation in the 
Sun (lRud igeill98 ^ an d in simulations (e.g.. lMiesch et alJl2006l : 
IWarnecke et alj|2013ah . 


3.2. Differential rotation and meridional circulation 


As shown in IWarnecke et al.l (1201 3al) . simulations with a coronal 
envelope are more capable of reproducing a spoke-like differen¬ 
tial rotation profile than runs without a coronal envelope using 
similar parameters (iKapvla et alJl2013l) . In this work, the param¬ 
eters of the runs are nearly identical, so we can isolate the influ¬ 
ence of the coronal envelope. In Fig. |4]we plot the angular veloc¬ 
ity, O = Qo ± sin 6, in the meridional plane and in Fig. |5]as 
a function of radius at two latitudes. In all of the runs the equator 
rotates faster than the poles, si milarly to the other simulations 
with s imilar Coriolis number s ([BrowrieLaL ^1 ^pvla et all 
1201 31 : IWarnecke et al.l l2013at lAugustson et al.ll2()15h . Further- 
more, all runs possess a local minimum in the angular velocity 
O at mid-latitud es, which has been shown to facilitate equa tor- 
ward migration (IWarnecke et al.l 1201 41: iKapvla et al.ll201^ . In 
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(r//2)sin0 (r//2)sin0 


Fig. 4. Angular velocity Q.(r,9)/Qo for Runs A1 and A3 {top row). 
Runs A3t and A5 {middle row) and Runs B1 and B3t {bottom row). 
Here, n = Oo + sin 6 is the local rotation rate. The dashed lines 
indicate the surface {r = R). 


runs without coronal envelopes the rotation profile has a simi¬ 
lar structure; the contours of constant angular velocity show a 
strong alignment with the rotation axis, following the Taylor- 
Proudman theorem. Run A1 possesses the strongest differential 
rotation. In particular at the equator the surface rotates faster than 
in the other runs. Runs Ale and Alc2 show a similar radial de¬ 
pendency, but differential rotation is weaker than in Run Al; see 
Fig- Ha) and (b). 

In the runs with a coronal envelope the Taylor-Proudman bal¬ 
ance is broken and the rotation profile has a more spoke-like 
shape; see Fig. |4] The differential rotation and the local mini¬ 
mum at mid-latitudes in these runs is much weaker. The min¬ 
imum also occurs at a greater depth than in the runs without 
a corona. The dependence on the size of the coronal extent on 
differential rotation is weak - similarly to what is seen for the 






r/R 


Fig. 5. Differential rotation D(r, d)/Do as a function of radius for all 
runs of Set A at d = 75° (15° latitude) (a), and at the equator 6 = Tijl 
(b), and Set B at mid-latitudes 6 = 75° (c), and at the equator 6 = njl 
(d). The inlay in (b) shows the angular velocity near the surface for 
Runs A2, A3, A3t, A4, and A5. 


latitudinal temperature distribution; see Fig. [3 The cooling pro¬ 
file has a minor influence on the rotation profiles such that in 
Run A3t the contours of constant rotation are slightly more ra¬ 
dial than in Run A3 In Run A5, the differential rotation is even 
weaker with more radial contours of rotation. 

The overall rotation profiles of Set B are similar to those of 
Set A and also show a local minimum of Q. at mid-latitudes al¬ 
though it is weaker and at a greater depth; see Fig. |4] Also here, 
the presence of a coronal envelope leads to more radial contours 
of rotation. However, the amount of differential rotation is re¬ 
duced only for Run B3t, but not for Run B3. 

If we compare the d ifferential rotation profiles with those of 
IWarnecke et alJ (l2013ah . we see that the contours of constant ro¬ 
tation are more cylindrical here. In particular, moderate rotation 
runs had more spoke-like differential rotation than rapidly ro¬ 
tating ones which is not seen in the present work. Therefore, 
we relate the more spoke-like rotation profiles in IWarnecke et al.l 
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Fig. 6. Meridional circulation in terms of the mass flux pu„, for 
Runs Al and A3t. The dashed lines indicate the surface (r = R) and 
the red solid line the inner tangent cylinder. 


(12013ah to a lower density stratification as the other parameters 
are similar. 


For most of the runs (A2, A3, A3t, A4, Bl, B3, and B3t) 
the maximum of rotation at the equator is actually below the 
surface, indicating a near-surface shear layer with negative ra¬ 
dial shear. The logarithmic gradient of rotation, dlnQ/dlnr, is 
around -0.2 near the surface for Run A4, and -0.15 for Runs A2, 
A3, and A3t. This is much weaker than th e value for the Sun, 
which is dlnQ/dlnr -1 for all latitudes (iBarekat et alJl20l4 
l2016l) . In all runs of Set B, except for B Ic, the near-surface shear 
region is more extended than the ones in Set A, and the gradi¬ 
ent is stronger; dlnQ/dlnr reaches values of -0.8 for Run B3 
and -0.5 for Run B3t. It is expected that for runs with lower 
rotation rate, a near-surface shear layer is stronger due to the 
weaker influence of the Coriolis force near the surface; see Sec¬ 
tions 13.31 an d 13.41 In agreement with A effect theory (iRiidigeil 
Il980l Il9^ . the double-logarithmic gradient should only be 
close to -1 v ery near the surface where the l ocal Coriolis num¬ 
ber is small (iKitchatinov & Rudigei 120051 : iKitchatin^ l20 1 31 : 
iRudiger et alJl201^ Kitchatinovl201f ). While this is true for the 


Sun, it is not the case in our simulations owing to limited strati¬ 
fication. 


To investigate the influence of the cooling profile on the tem¬ 
perature variation and differential rotation, we have performed 
two additional runs, in which we either increased or decreased 
the cooling luminosity compared to Run Alc2 (not shown). A 
decrease of the cooling luminosity by a factor of two leads to 
a shift of the temperature minimum at higher radii. The mean 
temperature increases slightly, resulting in a higher density at the 
surface and a decrease of the density stratification in the convec¬ 
tion zone. An increase of the cooling luminosity has the opposite 
effect, leading to a temperature minimum at a greater depth, and 
a lower temperature and density in the convection zone. Weaker 
cooling causes stronger differential rotation, especially at higher 
latitudes, while stronger cooling does not show a significant ef¬ 
fect. The gradient d In Q/d In r at the equator becomes more neg¬ 
ative with a weaker cooling. 

Differential rotation is also generated in the coronal en¬ 
velopes. Below r - 1.01 R, the differential rotation follows the 
rotational behavior of the convection zone. Above r = 1.01 R 
the plasma rotates nearly uniformly near the equator with a ro¬ 


tational speed close the Qq. The mid-latitudes rotate faster than 
the equator and at high latitudes the coronal envelopes decrease 
to slower rotation. All runs ex hibit strong cylindrical and ra¬ 
dial shear. This is consistent with lWarnecke et ^ (l2013ah . where 
runs with lower stratification showed a similar behavior. The 
change of behavior from below to above r - 1.01 R is caused by 
the change in temperature and density stratifications and not by 
the viscosity profiles. This can be seen by comparing Runs A3t 
and A5, where the same temperature and density profiles are 
used, but a much lower viscosity is applied in Run A5 than in 
Run A3t. 

In Fig. |6] we plot the meridional circulation in terms of the 
mass flux pu,„ in the meridional plane, where = {ur,U8,0) 
is the meridional flow. The meridional circulation has a multi¬ 
cellular structure in all runs. Near the equator at the surface the 
flow is poleward, but it can become equatorward at high lati¬ 
tudes; see Fig.|6] The strongest contribution to the mass flux car¬ 
ried by the meridional circulation occurs within the bulk of the 
convection zone. There the flow is aligned with the rotation axis 
and streaming toward the equator along the inner tangent cylin¬ 
der and toward higher latitudes further away from the rotation 
axis. These mass flows seem to stream toward the local minima 
of Q at mid-latitudes. From there, most of the runs develop a 
flow toward the equator following the 6 direction. The stronger 
meridional flows in Run Al are due to the higher density, see 
Fig. [TJb), while the actual flow is quite similar in all runs of 
Set A, see Fig.|8lh). The runs of Set B generate stronger merid¬ 
ional ci rculation, similarly to what was found in IWarnecke et al.l 
(l2013ah . At these rotation rates, slower rotation leads to an in- 
crease of mer i dional circulati on as found in mean-field models 
(Kohlei 19701: Rhdigeii 19 89h and nume rical simulations (e.g., 
iBrown et al .112008 : Augustson et al.ll^012h . In general the merid¬ 
ional flow pattern does not change due to the influence of the 
coronal envelope. 


3.3. Reynolds stresses and A effect 

Differential rotation and meridional circulation in the Sun and 
other stars is generat ed by the inter action of turbulent con¬ 
vection and rotation (iRiidigeii Il989t) . Reynolds stresses be¬ 
come anisotropic due to an angle between the gravity and 
the axis of rotation. The non-diffusive contribution of the 
Reynolds stress tensor Qij can be expressed via the A effect 
which produces equatorial acceleration if the a ngular momen¬ 
tum transport is directed equatorw ard (see e.g., lRudigeilll989l : 
IKitchatinov & Rudigeiill995Ll2005h . It has been recently shown 
that there is strong evidence for the A effect operat ing in the Sun, 
causi ng the observed rapidly rotating equator (lRudiger_et^ 
2014) and the latitude-independent surface shear (IKitchatinov 


20161 ). 


We calculate the three off-diagonal components of the 
Reynolds stress tensor Qyg - u'.Ug, Qr^ = Qg^ = 

Both Qr^, and Qg^ contribute to the angular momentum balance 
and their non-diffusive parts are associated with the horizontal 
and vertical A effects, respectively. Even though Qrs does not 
directly contribute to angular momentum transport, it has been 
argued to be important in gener ating a near-surfac e shear layer in 
global convection simulations (iHotta et al. The Reynolds 

stress es can be written as (e.g., iRiidigej 119891: iRiidiger et al.l 
120131) 


dQ 

= Ay smOD.- y^r sin 6 —, 
or 


(14) 
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Fig. 7. From left to right; Off-diagonal components of the Reynolds stress Qrg, Qr^p and Qg^ normalized by v,fio (first three columns), the turbulent 
viscosity in terms of molecular viscosity v,/v (fourth column) and the three components of the A effect Am, Ay, and Ah (fifth to seventh column) 
normalized by v, for Runs A1 (top row), A3t (middle) and A5 (bottom). 


= AhCos 0Q-vtsine-^, (15) 

where Ay and Ah are the vertical and horizontal components of 
the A effe ct and Vt is the turbul ent viscosity (assumed isotropic). 
Following lKanvla et al.l (l2014h we approximate Vt as 

A = ^u'nasO^MhjHp, (16) 


where aMLT = 5/3 has been assumed for the mixing length 
parameter, and Hp - -{d\Tip{r,0)ldrY^ is the pressure scale 
height. 

A meridional A effect also exists jPufkkinen,et^ 1^93 ’ 


Rieutord et al.lll994HKapvla et al.ll^04t iKapvla & Brandenbur 

20081 ) which is related to Qyg via 


In the earlier definitions, no Q. factor was included in the first 
term because, on theoretical grounds, one expects the rotational 
effects on the meridional part of the Reynolds stress Qrfl to b e 
proportional to sin 6 cos 6 (lRudigeilll989t iRieutord et ^.lil994l) . 
However, to obtain the same units for the coefficient Am as for 
the other components of A, we include here an Q factor. In Ap¬ 
pendix [A] we give a simplistic derivation for these coefficients, 
which confirms that Qre is not only proportional to a higher 
power of Q than the other two components of A, but that it 
also picks up a contribution proportional to Heir. However, under 
nearly isotropic conditions, the contribution of Am should have 
the same sign as ug', see Appendix lAl In our simulations this is 
indeed the case; see Fig. |7] 

We can now solve Equations (fT5l l- (fT7l i for Am, Ay and Ah 


Qr 0 - Am sin 8 cos 6Q. — Vi 


1 dur d 



(17) 


Am - 


Qre Vt 

sin0cos0Q sin0cos0Q 


1 dur d /ue\ 

7 "^ 


( 18 ) 
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Fig. 8. Off-diagonal components of the Reynolds stress Q^g (a), 2,.^ (b) and Qg^ (c) normalized by v,Oo. the turbulent viscosity in terms of 
molecular viscosity v,/v (d) and the three components of the A effect Am (e), Ay (f) and Ah (g) normalized by v, as well as the meridional flow 
“e/“rras (h) for Set A in the northern hemisphere at 15° latitude. The thin black lines indicate the zero value and the surface (j = R). 


Qi<t> 5 In Q 

-h Vtr —-—, 

sin0Q dr 

(19) 

Qg^, 5 In Q 

-h Vt tan 9 —-—. 

cos 6Q 89 

( 20 ) 


We plot the Reynolds stresses {Qre, Qrip, and Qg^), the tur¬ 
bulent viscosity Vt, and the three components of the A effect 
(Am, Av and Ah) in the meridional plane for Runs Al, A3t, 
and A5 in Fig. |2]and as a latitudinal cut (15° latitude) for all 
runs of Set A in Fig. [ 8 ] For Run Al, the Re ynolds stresses show 
the us ual behavior of rotating convection (iKanvla et al.ll20lil 
l2016cl) . Qrg and Qg^ are antisymmetric over the equator. In the 
northern hemisphere, Q^g is negative at the surface and posi¬ 
tive deeper down, but only outside the inner tangent cylinder. 
The latitudinal v ariatio n of Qrg agrees wit h that f ound both by 
iPulkkinen et al.l (Il993h and iRieutord et al.] (Il994t) . Qg^ is posi¬ 
tive in the northern hemisphere near the surface and weakly neg¬ 


ative deeper down at l ow latitudes. A s expected, Qr^ is symmet¬ 
ric over the equator (lRudigeiill98C)t) . This stress component is 
mostly positive (negative) at low (high) latitudes with a small 
negative region at the equator. The meridional structure of all 
R eynolds stress comp onen ts agrees qualitati vely with Run A 6 
of iKapvla et alJ (1201 Ih and iHotta et alJ (12015h with significantly 
lower and higher density stratifications, respectiv ely. However, 
the pe ak values are half of those in Run A 6 of iKapvla et al.l 
(1201 Ih . which is likely due to the slower rotation in their study, 
as the stresses are known to be que nched for faster rotation (e.g., 
lRudigeiiFl989l:lRudiger et al.ll2013h . 

The turbulent viscosity Vt has a maximum at the equator and 
towards the bottom of the convection zone. The value of Vt/v is 
close to Reynolds number Re = 35. Am is non-zero only at low 
latitudes, being positive near the surface in the topmost 5-10 
Mm, negative a bit deeper and mostly positive with a lower am¬ 
plitude even deeper down. Ay shows strong alignment with the 
rotation axis, being positive (negative) outside (inside) a cylin- 
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drical radius of 0.85 R, with a negative region near the surface 
at the equator. Ah has a similar structure, but shows a concen¬ 
tration near the equator above 0.85 R, where Ah is three times 
stronger than Ay. A simila r structure, bu t with three times lower 
values , has been found bv lKanvla et alJ (l2014ll and iKarak et alJ 
(12015h using the same technique, but for runs that rotate slower. 

The presence of a coronal envelope significantly alters the 
Reynolds stresses and therefore Vt and the A effect; see Figs. |7] 
and [8] The Reynolds stresses lose their alignment with the ro¬ 
tation axis. Qrg changes sign and is positive (negative) at lower 
la titudes in the north ern (southern) hemisphere, which is similar 
to lHotta et al.l (l2015h . This behavior has been found for all runs 
with an extended coronal envelope (A2, A3, A3t, A4, and A5), 
whereas in Runs Al, Ale, and Alc2, Qrg is negative close to the 
surface in the northern hemisphere; see Fig.[8la). This seems to 
confirm the presence of a correlation between a positive (neg¬ 
ative) value of Qrg in the northern (southern) hemisphere and 
the generation of near-surface negative shear. The inclusion of a 
coronal layer changes Qr^ such that the minimum at the equa¬ 
tor disappears and the overall magnitude of the stress is reduced 
by roughly a factor of two; see Fig. [T] At 15° latitude, all runs 
show a similar behavior with the main variation being the loca¬ 
tion of the maximum near the surface depending on the depth of 
the cooling layer; see Fig. Ob). Qg^ changes similarly as Qr^. 
The maxima near the surface are again shifted deeper in the runs 
with a deeper acting cooling layer; see Fig.Qc). The pro files of 
Qr^ and Qg^ of Run A3t are similar to lHott^t al.l (l2015h . 

The turbulent viscosity Vt varies less as a function latitude in 
runs with a coronal envelope; see Fig.l?] FigureOd) shows that 
Vt is reduced at the bottom of the convection zone and near the 
surface in the runs with a coronal envelope which is consistent 
with the lower turbulent velocities in those cases; see Fig.O^). 
In Run A5, the profile is similar to the other runs with coronal 
envelope, but Vt/v is higher due to a lower value of v. 

All of the A-coefficients are reduced in runs with a coro¬ 
nal envelope; see Figs. |7]andOe)-(g). This is consistent with a 
weaker differential rotation in these runs; see Fig. |5] Similarly 
as Qrg, also Am changes sign in the equatorial regions in the 
runs with a coronal envelope. This is consistent with the mean 
meridional flows changing sign in these runs as can be seen in 
Fig. [Sh). Ay does not change significantly apart from the re¬ 
duced amplitude. The vertical A effect is thought to be responsi¬ 
ble for generating radial shear. This is consistent with a stronger 
Ay producing a stronger radial differential rotation in Run Al 
and with a correspondingly weaker Ay and weaker radial differ¬ 
ential rotation in Runs A2, A3, A3t, A4, and A5; see Fig.|5] The 
horizontal A effect shows a more broad maximum as a function 
of latitude in the near surface layers in Run A3t in comparison 
to Run Al, whereas the amplitude is reduced by a factor of four. 
Figure [8]reveals that the values of Ah near the surface are lower 
in the runs with a coronal envelope or a cooling function pene¬ 
trating the surface, consistent with weaker latitudinal shear. The 
changes due to the coronal envelope are even more pronounced 
in the anisotropy parameters Am, Ay, Ah, which are related to 
the A effect; see Appendix iBl and Fig. IB.II Am and Ay change 
from being mostly positive to mostly negative (also within the 
coronal layer), whereas the changes in Ah are not as strong. The 
results for the more turbulent Run A5 are very similar with those 
of Run A3t; see Fig. |7] 

In Set B, where rotation is slower, the rotational influence on 
convection is reduced by roughly a factor of two. However, the 
only major difference to the runs in Set A occurs for Qrg, where 
the sign at low latitudes near the surface is positive already with¬ 
out a coronal envelope, although the difference in the meridional 





Fig. 9. The two dominant terms of Equation 1211 ) and their difference 
AL = rsmddSl^jdz + (VL x for Runs Al, A3t, A5, Bl, and B3t 

in the northern hemisphere at 15° latitude. The thin black lines indicate 
the zero value and the surface (r = R). 


A effect is rather minor. The main effect of the coronal enve¬ 
lope on the Reynolds stresses as well as the components of the 
A effect is their reduced amplitude similarly as in Set A. 


_ A ccording to theory (e.g., Ktehatino y & Riidig^ 1 19951 

120051 ; iRiidiger et al.l 120141 ; lKitchatinovll2016l) . the nei^surface 
shear layer in the Sun is caused by the vanishing horizontal A 
effect and the sole contribution of a negative vertical A effect. In 
our simulations, we do indeed And a weaker Ah in some of the 
runs, where d In Q/d In r is negative, see Figs.|7]and|3g), but for 
most of the runs the relation is inconclusive. In Section we 
will investigate this in more detail. 


In this section we have not discussed the influence of the 
magnetic held on dif ferential rotation g enerators. This has been 
discussed in detail in iKarak et al.l(l2015h for similar simulations, 
but without coronal envelope. They found that the large-scale 
and turbulent Maxwell stresses are around two orders of mag¬ 
nitude smaller than the Reynolds stresses and therefore do not 
influence the angular m omentum transport si gnificantly. How¬ 
ever, the recent study bv iKapvla et al.l (12016cl) has revealed that 
the turbulent Maxwell stresses in similar setups without coronal 
envelopes are of the same order of magnitude as the Reynolds 
stresses at high magnetic Reynolds numbers. The current runs 
are in the intermediate range in ReM where the total stress is al¬ 
ready affected but the qualitative character of the hydrodynamic 
results remains unchanged. 
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3.4. Thermal wind balance 


The Taylor-Proudman balance can be broken by a non-zero 
baroclinic term in the mean azimuthal vorticity equa¬ 
tion, also know as the thermal win d balance or meridional 
circulation evolution equation (e.g., Brandenburg et alJ 119921: 


iKitchatinov & RUdigedll9^IW^necke et al.ll2013ah. 


dtL>s dQ.^ 

—— - rsmO— -h 

dt dz 


|vr X Vs] 


V X - V ■ pu'u' 
P 


, ( 21 ) 


where djdz - cosOdjdr — sin0<9/(90 is the derivative along 
the rotation axis and ru = V x « is the mean vorticity. We ne¬ 
glect here the contribution of the Maxwell stress arising from the 
correlations of the fluctuating magnetic field and those from the 
mean magnetic held. The contributions of meridional fl o ws turn 
out to be small and can be neglected. IWarnecke et al.l (l2013ah 
could show that the baroclinic term, that is the second term in 
Equation (12111 . is responsible f or spoke-like differentia l rotation 
profile. Following the study of lWarnecke et al.l (l2013al) . we plot 
the two dominant terms on the right hand side of Equation (l2Tli 
together with the residual AL = rsin0(9f2^/(9z -i- |^VLx Vsj 
in Fig. 0for Runs Al, A3t, A5, Bl, and B3t. The meridional 
distributions of |^VL x and rsin05Q^/(9z of Runs Al and 

Runs A 3t are very similar to Figures 5 and 9 of lWarnecke et al.l 
(l2013ah . respectively. For most of the convection zone the two 
terms match well and AL is close to zero. In the upper 0.1 R, AL 
is non-zero. In the case of Run Al, AL is positive in the region 
above r - 0.94 R. However, right at the surface, it is negative, 
probably because of boundary effects. By contrast. Run Bl de¬ 
velops regions, where AL is positive (0.90 R < r < 0.94 R) and 
where AL is even more strongly negative (0.95 R < r < 0.99 R). 
In Runs A3t, A5, and B3t, the overall magnitude of terms in 
the thermal wind balance is lower, resulting in a significantly 
smaller residual AL. However, AL seems to be close to zero in 
all three runs and becomes non-zero only near the surface; see 
Fig. 0 Thus, there is no clear change in behavior depending on 
rotational influence or Prandtl numbers in the coronal runs. As 
the meridional circulation is stationary, apart from cycle depen¬ 
dent variations, the left-hand side of Equation (12111 is on average 
zero. Therefore, in the region close to the surface the contribu¬ 
tion of the Reynolds stress, that is the third term on the rhs of 
Equation (l2Tli . plays a more important role and will be related to 
the residual AL. We can rewrite this term 


AL 


V X ( - V pu'u' 


- [Vp X V ■ pl/'H'j 


= -:^[(9rp[V pH'H']^- -(9flp[v pH'H']^ 


( 22 ) 


The latitudinal component of the divergence of p u'u' is given by 


[V-pM'M']^ = —3r[t^pu'^u'g) 


pu'^g 


I , _, cot 6 pu'.u', 

—^dg{r sin Opu'gu'g) -(23) 

in « ' a nj 


r sin0 






Fig. 10. The three dominant contributions of the Reynolds stress Qrs 
llSl, Qee (126) and ^ as well as their sum for Runs Al, A3t, A5, 
Bl, and B3t in the northern hemisphere at 15° latitude. The thin black 
lines indicate the zero value and the surface (r = R). 


where u'.u'. — m'm' and 5/50 = 0 because of the azimuthal mean. 
This means that the contributions related to u'.u',, and u'„u', are 

T 0 ’ 0 (p 

zero. Furthermore, we And that terms related to u'.u', are too 
small to have a strong effect. In Equation (l23b the first three 
terms and in Equation (l24l) only the first term have a strong con¬ 
tribution to the thermal wind balance. We summarize them in the 
following three expressions: 


--5.P 

P 




2--TT\ 

r pn^M') + —— 




r9. 


where the first term on the Ihs is the dominant one. 


(25) 



and 


(26) 


and the radial component of the divergence is given by 

[V'Ph'h'I = \-drir~pu',.u',) + —^—5fl(rsin0pM'M' 

L Jr 1-2 \ ' rsm 6 ^ ■ 

pu'gu'g+pu'^u'^ 

r 


(24) 



(27) 


In Fig.[T0l we plot all three terms and their sum for the same runs 
as in Fig.0 In Run Al, their sum is positive, resulting in a neg¬ 
ative contribution of the Reynolds stresses to the thermal wind 
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Fig. 11. Inclination of the magnetic field near the surface plotted as 
a 2D histogram over radius r/R for Runs Al, Alc2, A3, and A3t. 0° 
means fully radial and 90° fully horizontal. The red line indicates the 
average for each radii. 


the field and the radial direction is distributed nearly uniformly, 
as indicated by a mean of the absolute angle being close to 60°. 
The field is fully radial just in the two uppermost grid points; see 
Fig- HU In Runs Ale and Alc2 the field distribution is the same, 
except that the magnetic boundary condition is placed higher. 
For the runs with coronal envelope, the field is less radial at the 
surface than in Run Al. Furthermore, the field is more horizontal 
in the coronal envelope than in the convection zone. The compar¬ 
ison of Runs A3 and A3t suggests that the cooling layer, which 
reaches deeper in Run A3, allows the field to become less radial 
at and above the surface. This leads us to conclude that the radial 
boundary condition does not significantly affect the structure and 
inclination below r - 0.98 R. 


IWarnecke et al.l (12016bh find turbulent radial downward 
pumping from a similar run, which can cause radial alignment 
of the magnetic field, but only close to the surface (r > 0.98 R) 
and at low latitudes, which is consistent with our findings. In 
local convection simulations, it has been found that turbulent 
pumping is a major ef fect causing horizontal magnetic field to 
move downward (e. g., Nordlund etmil992t iTobias et al.lll998t 


___ f e _ 

lOssendriiver et al.ll2002 ). It is also believed that this effect can 


cause the field near the surface to become dominantly radial 
(ICameron et al.ll2()TU) and is therefore the m ain reason for the 
succe ss of surface flux transport models (see iMackav & YeatesI 
1201 2 i for a review). 


balance; see Equation (l2Tli . In Run Bl, these terms are positive 
around r = 0.94 R and negative closer to the surface. Qrr is in all 
runs negative giving a positive contribution to the thermal wind 
balance. The runs with a coronal envelope seem to have a weaker 
and negative contribution from the Reynolds stresses. This is 
probably related to the low density stratification near the surface. 
If we compare the plots of Fig. [T0]with Fig. |3 there seems to be 
some agreement of the residual AL with the contributions of the 
Reynolds stresses Qrg, Qgg, and Qrr- However, there seems to be 
a shift in radius, which might be related to still missing contri¬ 
butions, which we could not identify in the present work. This 
might be related to the fact that in our simulations, in particular 
near the surf ace, V ■« is not zero as assumed in m any mean-field 
models (e.g.. lKitchatinov & Rudigei1ll995Ll2005h . Furthermore, 
we have here neglected the contribution of the Maxwell stresses, 
which might also expl ain this inconsisten cy. However, we can 
confirm the result from iHotta et al.l (12015h that Qrg is important 
near the surface and gives an important contribution to the ther¬ 
mal wind balance near the top boundary. However, this contribu¬ 
tion is not directly related to the stress Qrg itself but to its radial 
gradient. Furthermore, also Qgg and Qrr have a large contribu¬ 
tion to the thermal wind balance in terms of Qgg and Qrr- Similar 
to Qrg in Fig. 0 Qrg changes sign in the runs where we have 
identified a near-surface shear layer; see Section [3721 This effect 
is even stronger near the equator, where it could play a role in 
generating the near-surface shear layer. 

3-5- Radial field boundary condition 

In the case of Runs Al and Bl, we employ a radial field con¬ 
dition at the surface, r - R, whereas in the other cases the ra¬ 
dial condition is enforced either at r = 1.01 R (Runs Ale, Alc2 
and Bl) or high above the convection zone (Runs A2, A3, A3t, 
A4, A5, B3, B3t). This affects the magnetic field distribution 
and therefore the dynamo in the convection zone. In the major¬ 
ity of the upper convection zone of Run Al, the angle between 


3.6. Cyclic dynamo solutions with dynamo wave propagation 

In all runs of both sets, convective motions, overall rotation, and 
their interaction contribute to generating a large-scale magnetic 
field. Most of the runs produce cyclic magnetic fields in the 
saturated phase. In Fig. [12] we plot the time evolution of the 
mean toroidal magnetic field for a selection of runs. To in¬ 
vestigate the cause of the propagation di rection of the magnetic 
fiel d, we apply the same technique as in IWarnecke et al.l (l2014tl 
and iKanvla et al.l (12016ah and we calculate the pr opagation di¬ 
rection using the so-called Parker-Yoshimura rule (iParked! 19551 : 
lYoshimu^l 1975h . 

As discussed in IWarnecke et al.l (l2014ll in detail, Runs Al 
and Alc2 show a solar-like equatorward migration of the mean 
field. This is probably caused by an equatorward propagating dy¬ 
namo wave generated by the region of negative shear in the mid¬ 
dle of convection zone; see Fig.|4|and bottom row of Fig.[T2| We 
have confirmed this interpretation using all tur bulent transport 
coeffic ients obtained with the test-field method (IWarnecke et al.l 
12016bl) . To reiterate. Runs Al and Alc2 have two dominant dy¬ 
namo modes; equatorward and poleward migrating branches at 
mid-latitudes and high latitudes, respectively; see top row of 
Fig. IT 2 I The dynamo wave seems to be formed in the middle of 
the convection zone and propagates toward the bottom and top 
of the convection zone; s ee bottom row of Fig. IT2] The cycle pe¬ 
riod is around five years (IWarnecke et al.l2014li . Near the surface 
there is a weaker mode, which has a cycle period of 1.5 years 
that propagates poleward at mid-latitudes and vanishes at higher 
latitudes. At the bottom of the convection zone exists another 
dynamo mode, see Fig.[T2| which might cause a long-term vari¬ 
ability. Such competing dyna mo modes have been analyzed in 
detail in iKapvla et al.l (l2016ah . Furthermore, the magnetic field 
in Run Ale shows no equatorward migration. As described by 
IWarnecke et al.l (l2014l) . this is most likely caused by the much 
weaker negative radial gradient of Q. in the middle of convection 
zone, also visible in Fig. [5] a) and the suppression of turbulent 
convective motion near the surface, as discussed in Section ITTI 
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Fig. 12. Mean toroidal magnetic field evolution as time-latitude (butterfly) diagram, plotted at a radius r = 0.98 R {top row) and time-radius 
diagram plotted at a 25° latitude (middle row) in kG during a 20 yr interval in the saturated stage for Runs Al, Alpc, A3, A3t, and Bl. For Run Bl, 
we plot the butterfly diagram at radius r - 0.84 R and its time interval is 40 yr. The dashed horizontal lines mark the equator (9 = nil) and the 
radii r = R, r - 0.98 R and r - 0.85 R, respectively. Propagation of the mean magnetic field for the same runs (bottom row). Color coded is 

S lotted during the saturated sta ge to gether with white arrow s showing the direction of migration ^mig(F S) = -ae^ x VQ of an afl dynamo wave 
Parke jfT95.5riYoshimurall 1 97.5h : see lWamecke et al](l2014h . We suppress the arrows above r = 1.05. The black solid lines indicate isocontours of 
at 2.0 kG. The dashed white lines indicate the surface (r = R). 


In this paper, we focus on the differences between the runs with 
and without a coronal envelope. 

Before we discuss the runs with a coronal envelope, we in¬ 
vestigate how the magnetic evolution changes, if we change the 
radial field condition to a perfect conductor condition at sur¬ 
face (Run Alpc). This is sufficient to cause the dynamo modes 
to change. Close to the surface the toroidal magnetic field be¬ 
comes the strongest, reaching values of more than 10 kG which 
are more than two times higher than the maximum values in 
Runs Al and Alc2; see top and middle row of Fig. [121 There the 
predicted migration direction is poleward, which agrees with the 
actual migration near the surface; see bottom row of Fig.[T2l Fur¬ 
thermore, we find indication of a poleward migrating subdomi¬ 
nant dynamo mode similar as in Runs Al, Ale, and Alc2. How¬ 
ever, we find an equatorward propagating mode in the middle of 
the convection zone, which is most likely related to the smaller 

and weaker concentration of found in the middle of the 

0 

convection zone which is predicted to propagate equatorward. 
It seems that the hydro-thermal setup tends to produce strong 
toroidal magnetic fields near the surface, but the radial field con¬ 
dition at the top boundary of Run Al, Ale, and Alc2 prevents 
this. The toroidal magnetic field is concentrated at mid-latitudes 
and produces nearly no field above 35° latitude. This is very sur¬ 
prising, given that the other runs without a coronal envelope pro¬ 
duce strong magnetic fields at high latitudes. This absence of po¬ 
lar field in Run Alpc suggests a relation between the polar field 
and the radial field at the surface. The relation might be interest¬ 


ing in view of the Babcock-Leighton dynamo framework (e. g., 

lBabcocklll961tlLeightonll964t iDiknati & Charbonneaul[l99^ . 

Adding a coronal envelope on top of the convection zone 
changes the magnetic field evolution significantly; see Fig. [T2] 
For the runs of Set A, the mean toroidal magnetic field migrates 
mostly poleward in a region between the equator and +40° lat¬ 
itude. At high latitudes the magnetic field is weaker, but shows 
a tendency of equatorward migration, in particular in Run A3t 
this equatorward propagation is clearly visible. In the middle of 
convection zone, the equatorward propagation also dominates at 
lower latitudes. As a main difference compared to Run Al, the 
mean toroidal field occurs close to the surface instead of in bulk 
of the convection zone. At this location, the radial shear is posi¬ 
tive causing a dynamo wave to propagate poleward, see bottom 
row of Fig. [12] This agrees with the magnetic field propagation. 
Furthermore, we find some indication of predicted equatorward 
migration in the middle of the convection zone, however this is 
not as clear as in Run Al. We do not find any clear difference 
between Runs A2, A3, and A4 (only Run A3 is shown). This im¬ 
plies that the size of the coronal envelope is not important for the 
magnetic field evolution. However, in Run A3t with a modified 
cooling profile the field is weaker in the band near the equator 
and more concentrated at or even above the surface in compari¬ 
son with the other runs with coronal envelopes. Furthermore, the 
equatorward migration near the poles is more pronounced. 

Also the migration periods of the large-scale fields change 
due to the coronal envelope. The cycle period of the poleward 
migrating field in Runs A2, A3, and A4 are around two years. 
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which is sh orter than in Runs A1 and Alc2, where it is around 
five years (IWarnecke et al.l l2014li . However, the equatorward 
branch near the poles seems to appear only every second pole- 
ward cycle in Runs A2, A3, and A4, which give a similar period 
as in Runs A1 and Alc2. In Run A3t, the magnetic field near the 
equator does not show a regular behavior with a clear cycle. The 
equatorward branch near the poles has a period of around two 
years. The magnetic field evolution of Run A5 shows the same 
features as A3t; poleward migration or quasi-stationary behavior 
at low latitudes and equatorward migration at high latitudes. 

In Set B with slower rotation, the magnetic field evolution 
shows a similar dependence on the cooling profile and the coro¬ 
nal envelope as in Set A. In Run B1 the magnetic field shows 
equatorward migration, in particular in the middle of the con¬ 
vection zone, as shown in Fig. The radius-time diagram is 
similar to that of Run Al, but due to the slower rotation, the cy¬ 
cle period is extended to around ten years. This is consistent with 
the predicted dynamo wave scenario; see bottom row of Fig. 
Analogously to Run Ale, the cooling layer in Run Blc causes 
the magnetic field to lose its equatorward migration mode. The 
other runs have mostly a stationary mode compared to their 
rapidly rotating counterparts, in particular at low latitudes. As 
in Set A, all migration directions of the magnetic field in runs of 
Set B agree with the predicted propagation. 

In general it seems as if the toroidal magnetic field tends to 
be strong near the surface in the runs with an extended coronal 
envelope (Runs Alpc, A2, A3, A3t, A4, B3, and B3t), but in 
the runs with the vertical field boundary condition (Runs Al, 
Ale, Alc2, Bl, and Blc) the boundary condition does not allow 
this. Overall, there is good agreement between the predicted and 
actual direction of migration of the mean toroidal field. 

The magne t ic field evolution is similar to that found in 
IWarnecke et al.l (l2013ah whose runs have more than two times 
higher SGS Prandtl numbers and lower stratification (p/psurf = 
14), but the other parameters (Re, PrM, Co) are comparable with 
the runs of t his work. In con t rast to this work in the rapid ro¬ 
tating runs of lWarnecke et al.l (12013ah the field near the equator 
seems to continue to propagate equatorward. In their slower ro¬ 
tating case the field becomes quasi-stationary in the saturated 
stage also at high latitudes. At this point it is not clear if these 
differences are related to the larger SGS Prandtl numbers or to 
the weaker stratification. 


4. Conclusion 

In this work we have studied the influence of the upper boundary 
on convectively driven dynamos in spherical wedges. For this 
purpose we have added a convectively stable coronal envelope 
of different sizes on top of the convectively unstable dynamo 
region. This coronal envelope effectively corresponds to having 
a free boundary as opposed to a stress-free radial field boundary 
condition used in man y earlier one-layer dyna mo simulations. 
We confirm the result of lWarnecke et akl (1201 3al) that the coronal 
enyelope has an influence on the dynamo region leading to a 
change in differential rotation and magnetic field eyolution. If the 
radius Rc of the coronal enyelope is just 1 % of the solar radius R, 
its influence is small and can entirely be related to small changes 
in the radial density and temperature profiles. If the size of the 
corona is larger {Rc >1.2 R), the influence is stronger. Howeyer, 
runs with coronal sizes extending higher than Rc = 1.2 R are 
nearly identical. 

Regarding the hydrothermal properties, the influence of the 
corona can be summarized as follows: (i) The radial mass flux 


across the free surface does not show any major difference be¬ 
tween runs with a small and an extended coronal enyelope. The 
radial mass flux in our simulations is too small to haye an influ¬ 
ence on the flow properties inside the conyection zone, (ii) The 
latitudinal temperature yariations due to rotation are significantly 
weakened due to the presence of the coronal enyelope. (iii) This 
seems to cause the differential rotation profile to become more 
spoke-like and weaker in runs with a coronal enyelope. (iy) This 
effect can also be seen in the change of the off-diagonal Reynolds 
stress components due to the coronal enyelope which can be ex¬ 
plained by a change in the A effect and the anisotropy. 

Furthermore, in the cases with G = 5 and a coronal enyelope, 
as well as in all runs with Q = 3, we find the generation of a 
weak near-surface shear layer at low latitudes. We haye related 
this generation to a change of sign in the meridional Reynolds 
stress tensor component Qrs near the surface in these runs. This 
component contributes to the meridional A effect and turns out 
to be non-zero in all simulations. Additionally, we haye shown 
that the radial gradients of Qrg and Qrr as well as the latitudinal 
gradient of Qgg are important near the surface to balance against 
differential rotation and the baroclinic term. A change of sign 
in these terms can be associated with the generation of a near¬ 
surface shear layer. 

The coronal enyelope seryes as a free top boundary for the 
magnetic field. We find that the dynamo properties are gener¬ 
ally strongly influenced by the choice of boundary conditions. 
This can be particularly important when the magnetic field is 
strong near the boundary, as was seen recently in connection 
with the l atitudinal bounda ry condition for simple mean-field 
dynamos (ICole et al.ll2016l) in the wedge geometry that we con¬ 
sider here too. In all of our simulations a toroidal field tends to 
form preferentially near the surface of the conyection zone but it 
is pushed down by a radial field boundary condition. Howeyer, 
with the radial field boundary condition, the field is only radial 
oyer a few grid points below the boundary; otherwise it is dis¬ 
tributed for all simulations i sotrop ically within the conyection 
zone. As in IWarnecke et al.l (l2014l) . we compare the migration 
of the mean toroidal magnetic field with the predicted propa¬ 
gation direction of the a£2 dynamo waye following the Parker- 
Yoshimura rule. It turns out that this rule can explain all the dif¬ 
ferent migration directions found in our simulations. This is a 
remarkable result giyen the yariety of the simulations. 

We must emphasize that the combination of a dynamo re¬ 
gion with a coronal enyelope is still far from realistic. This is 
mostly because of the rather low density contrast and the strong 
yiscous coupling between the two layers. We must therefore be 
aware of the possibility of artifacts, for example the occurrence 
of a differentially rotating coronal enyelope at higher latitudes 
may be an example. This is not normally expected to be the case 
in the Sun (iTimothy et al.lll97^ . although the obseryations haye 
only been limited to regions where one sees a rigidly rotating 
coronal magnetic fiel d. Indeed, recent work on coronal holes by 
iLionello et al.l (l2005l) claims that the rigid rotation is only appar¬ 
ent. 

Our work may haye an impact on our understanding of the 
properties of the conyection zone and the solar dynamo, be¬ 
cause changing the properties of the boundary and studying its 
influence teaches us about the physical properties and dynami¬ 
cal effects within the solar conyection zone. Howeyer, we also 
haye to keep in mind that our model, simila rly to other mod¬ 
els of spherical c onyection for the Sun (e.g., iBrun et akl l2004t 
iHotta et akllloBh . uses an SGS model for the unresolyed con- 
yectiye heat flux. In addition, these models represent the strong 
radiatiye cooling at the surface just by a relaxation term, which 
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ignores the effects of s trong driving of convection by w hat is 
known as entropy rain (ISprui ^ 1199'^ lBrandenburg||2016l) . This 
can have dr amatic effects on the ty pical length scales of stellar 


n the ty i 
^ 120161 ) . 


which in turn would affect 


convection (ICossette & Ras 
the differential rotation and dynamo properties. 

In our simulations, we have confirmed the influence of the 
latitudinal temperature distribution o n the differentia l rotation as 
described by mean-field m odels (e.g., Rudigel^l989^ and backed 
up by simulations (e.g., iMiesch et al.l l200« Warnecke et akl 
l2013ah . Only a small temperature difference between pole and 
equator is needed to cause the differential rotation to become 
more spoke-like. However, if the rotation causes too large lati¬ 
tudinal temperature variations, the differential rotation becomes 
dominated by the Taylor-Proudman theorem and has cylindrical 
contours. On the other hand, the strength of the differential ro¬ 
tation benefits from large latitudinal temperat ure variations. Fur- 
thermore, our work confirms the results of iHotta et al.1 (1201 5h 
in that Qre may be important for the generation of the near¬ 
surface shear layer. However, we went a step further and identi¬ 
fied the terms in the thermal wind balance which are important 
ne ar the surface. These results appear to be in conflict with those 
of iKitchatinov & Riidigetl (ll995Ll2()()5h . which explain the near¬ 
surface shear layer solely by the vanishing horizontal A effect 
near the surface. We confirm in some simulations that Ah goes to 
zero near the surface, but there we find no clear relation between 
Ah and a negative radial gradient of rotation near the surface. 
The associated Reynolds stress component also does not play a 
role in the thermal wind balance. Future helioseismic measure¬ 
ments may yield information about the presence and distribution 
of Qre in the Sun. 

The fact that the Parker-Yoshimura rule can explain the mi¬ 
gration direction of mean toroidal magnetic field found in our 
simulations has an impact on the interpretation of equatorward 
migration of the magnetic field in the Sun. Applying this rule 
to the Sun will lead to the generation of toroidal field in the 
near-surface shear layer, where negative shear can cause equator- 
ward migration. This would also imply that sunspots are formed 
near the s urface by a local flux conc e ntration mechanism as pro ¬ 
posed by Brandenburget^ J201 Ih. IStein & Nordlund j|2012h . 

I Warnecke et al.l ( 201 3bl 2016a ). and lKanvla et al.l ( 201 6bl) .l^- 
thermore, we find that in our simulations turbulent pumping is 
not strong enough to cause a preferred radial orientation of mag¬ 
netic field near the surface as in the Sun. This is possibly due to 
insufficient stratification in our simulations. However, a detailed 
analysis of the effect of tu rbulent pumping cruc i ally depends on 
the tes t-field method le.g. lSchrinner et al.ll2007t fWarnecke et al.1 
1201 6hh . 
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turbulent correlation time t, that is Rij - -Qijlr. Inserting this 
into Equation (IA.51 l and assuming a steady state, that is Qij - 
0 , using that the background velocity correlation is of the form 


u'.u'. - diag(M'^, u'^, u'^), we have 

Qrg = 2t(m'^ - Ugjr + ... , 

(A.6) 

Qrj, = Q sin 0 -h ..., 

(A.7) 


Qgj, = 2t(m^^ - u'g^) Q cos 0 H-..., 

(A.8) 


where has been replaced by Qrsinfl and gradient terms 
of Mg/r are assumed to vanish, so that Equation (IA.3II yields 
ii'g - -2{u'i.Ug+Uru'g)lr+... . We note that, while Qr^ and Qg^ are 
proportional to Q, the meridional component Qrg is proportional 
to Ug (= ug). If we were to allow ugfr to have non-vanishing gra¬ 
dients, we would have Qrg - t{2u'^- u'l^) Uglr+.... In that case, 

under isotropic conditions (m'^ = Ay = Ah = 0, but 

Am is non-vanishing and would be positive. Therefore, the non- 
diffusive contribution to Qrg would be negative for a poleward 
flow. This is in agreement with the profiles of Am and Qrg in the 
surface regions of our simulations; see Eig.|2] 

Appendix B: The anisotropy parameter 

The A effect is related to anisotropy parameters (lRudigedll980l) . 
defined as 


Appendix A: Meridionai A effect 


IM 


To obtain an expr ession for the A effect, we use the minimal 
tau approxi mation ( Blackman & Eield|2002, l2003h : see also the 
derivation in iKaovla & Brandenburg 12008 ). We denote partial 
time derivatives by a dot and compute 


Qii - ii'.u'. + u'.ii 

^ij I J I 


I J 


(A.l) 


in a non-rotating frame of reference (U - u+e^Q.Qr sin 9), using 
ii'r -2{ugUg + u'^U^)lr + ..., (A.2) 


/2 


/2 


.,'2 


,,2 


+ U'r^, 


+ U'r^, 



+ U'g^. 


(B.l) 

(B.2) 

(B.3) 


ii'g - -(u'^Ug + Uru'g)lr + 2u'^U^ cot6>/r + ..., (A.3) 

ii'^ = -2(M't/0 + Uru'^)lr-2(ugU^ + Ugu'^)cot0lr + ..., (A.4) 

where the three dots denote nonlinear and gradient terms of Ur, 
Ug, and M^/r sin 6 that will be neglected. Below we also consider 
the case where gradient terms of ug will be included. In other 
words, the base state corresponds to rigid rotation with a merid¬ 
ional flow Ug oc r sin 6, and gradients around this state are ne¬ 
glected. 

Inserting Equations (IA.2l) - (IA.4b into Equation (lA.ll) . we ob¬ 
tain an expression of the form 

Qij = UjklTk + Rij, (A.5) 

where Lijk is a rank 3 tensor (related to the coefficients of the 
A effect) and Rij stands for all remaining terms, in particular 
those that stem from the triple correlations. In the minimal tau 
approximation we replace those by a relaxation term with the 
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Fig. B.l. Anisotropy parameters Am, Ay, Ah for Runs Al {top row) and 
A3t {bottom). 
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